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ABSTRACT 


The interaction of land 5 kev copper atom (primaries) impacting the 
011 surface of copper crystals was studied using an "n-body type" 
computer simulation program. Each primary energy was studied using 
both the Bohr and Gibson II copper on copper potentials. Primary pene- 
tration and lateral drift or spread were investigated as a function of 
impact point. 

The n-body program has greater penetration values than the binary 
model used by Robinson and oon . Drift or spread (defined as wander 
if energy loss rate is low) occurs in the impact region of transition 
from hard interaction (near a lattice atom) to soft interaction (open or 
channel regions). Penetration and spread contours are presented. 

The guidance and assistance given by Associate Professor Don E. 
Harrison, Jr. of the U. S. Naval Postgraduate School in this study is 


greatfully acknowledged. 


ii 





TABLE OF CONTENTS 


Section Page 
ee lateoduction (7.4.5 Oy Oe Bane 2 ease ee eee 1 
2. Computer SimulationeStiudiles.34....288: 45 .60..o eee 6 
Spe ODJECtIVeS: "age. EERE RIS - Guin oul) cia ere eee 10 
A, Simulations Models . «ccc crewnee a 6 oa sles 6 a ces oi eee tl 
5. Procedure er re 7 ee ee aiRige aerate errors | 
ie Results ande@emelucloncwwee. ..- sr sect es wees ee 24 
7. Biblidgraphv.. a.) Ucn te aa eIeRRaIRE -si-0. 10 eames sa) 0 5s pommel 
Appendix 
[, Program Listing se. ar ote Renn Re neMer ete tetoner erste) cnet ner oe 
II. Program Discussion.......... Penne) oteiie tenets sl a 61 
III. Program Operating Instructions -----.-+ee-eeeee. coeeeesee 69 


Lod 





LIST OF ILLUSTRATIONS 
Figure Page 
Ly Integral Penetration Curves... oc descr ee eee ae 30 


2. Robinson and Oen "...(110_ channel trajectories 


onto (011) surfacesfor ivever Cu” ae we eee 36 
S (1T0) Face of Chanh2 Type Crystal... ae. ene ed 
4, (110) Face of Chana 1lype Unit Crystals... oer 38 
S. N-body Model Trajectories ....... PT eee Sse eee oes eee 
c. IMmtegral PenetationeGumvece errr ee eee 40 
7. Depth of Penetration and Maximum Scatter 

Gregg Reem Ch ee ee ee eS 21 
8. Impact Triangle, 14 and 12 point grid systems.......... 42 


9. Contour Curves. Gibson II potential and 
SUOQeEV primaries eres ce lc et ores sede Oe ee bac here 43 


10. Contour Curves. Gibson II potential and 
A OG a MAAS Sa oa Gace nics esses, 6 So alee: 6 0 SURE e oe “ale 44 


1). Contour Curves. Bohr potential and 5000ev primaries.... 45 


2x Contour Curves. Bu..r potential and 1000ev primaries .... 46 


3" Potential FunctionsG@urves....cssscse ce 4 atta levacs) coee een eee 47 
14, Impact Triangle, General ResultS .......c.ccocevscvccece 48 
Poy aw lGiosonslLPotential.5000ev, Primaries... snuwisdssansucdeweeG 





1, INTRODUCTION. 

The subject of crystal bombardment by energetic ions or atoms may 
be discussed from either external effects or internal interactions. 

Neither phenomenon is well understood, nor have adequate theories been 
devised to explain these interactions. One external effect is described 
by the sputtering ratio, which is defined as the ratio of ejected atoms 
from a crystal per impact atom. A more general term for crystal modi- 
fication caused by incoming energetic atoms is radiation damage. 

When incoming atoms or ions (called primaries) penetrate the crystal, 
a series of collisions may take place in which the primary energy is 
distributed to the crystal atoms. Often crystal atoms leave their sites 
and become moving atoms called secondaries. The orientation of the 
crystal, relative to the primary beam direction, is a critical parameter 
in any description of the interactions that take place. Deeper penetrations 
occur when primaries enter more open or transparent channels in crys- 
tals. 

The primary energy effectively categorizes the type of interaction 
which will occur. Heavy incoming atoms with relatively low energies 
usually interact with the atoms of the crystal through elastic collisions. 
Such atoms deviate markedly from straight line trajectories. High energy 


primaries interact by ionization and electron excitation of the crystal 


ae 








and, asa result, lose little energy. These atoms usually follow 
relatively straight line trajectories. te 

This study is concerned with the ranges and spread of lower ener- 
gy primaries as they dissipate energy through internal elastic inter- 
actions in the crystals. Range is definea as the depth of penetra- 
tion of the primary from the surface. For purnoses of this study total 

path distance and vector range from impact point to a specified point 
are not needed. Spread is defined as the radial displacement of the 
primary from the axis of its original impact line. 

These introductory sections provide a background in the recent 
experimental and machine calculation studies of low kev atoms and 
ions in crystals. The topics of sputtering and basic radiation dam~ 
age are not primary objectives of this thesis. However, several 
good sources of information are recommended: a basic text for intro- 
ductory study, by Dienes and Vineyard , 2 and more advance articles 


4,9 and Harrison, et. aie 


by Wehner, et. al. e Nelson, et. al., 
Although recent experimental range studies have been performed 
with heavy ions bombarding single crystal targets; the nature of these 
experiments and the lack of knowledge concerning potential functions, 
make the comparison between laboratory and simulation data diffucult. 


Until recently, the only theoretical potential functions available were 


* All footnotes refer to the Bibliography 





those from homonuclear atomic interaction. Thus, the computer pro- 
grams have used primaries and crystals of the same type, while in 
laboratory experiments the primaries and crystals have been differ- 
ent substances. 

Virtual leaders in the field of experimental range studies are J. A. 
Davies and his co-workers at Chalk River, Ontario. Their technique 
is both sensitive and reliable, but requires that radioactive gas ions 
be used as primaries. The approach entails: 


.-- bombarding an optically flat aluminum target with a mono- 
energetic beam of radioactive ions, and then measuring the 

depth of penetration of the incident beam by dissolving succes- 
ive uniform layers of aluminum from the target surface and measur- 
ing the amount of radioactivity in each layer. The ion bombard- 
ments were carried out using an electrostatic accelerator spec- 
ially designed for the purpose. The technique... (dissolving ex- 
tremely thin layers of known thickness)...consists essentially 

of two steps:electrochemical oxidation at constant voltage in 
aqueous ammonium citrate, followed by chemical removal of the 
anodic oxide film in a phosphoric acid chromium trioxide solution. 
Due to the highly protective nature of theanodic oxide film, this 
process permits highly reproducible surface layers of metal as 

thin as 37k to be removed. ...The process may be repeated as of-~ 
ten as desired, so that it is possible to obtain the complete range- 
distribution curve froma single bombardment. ’ 


Some of the earlier polycrystaline range studies performed by Davies, 


et. al. included .7 to 60 kev Na2 ions in aluminum,’ and 2 to 600 kev 


Kr89 ions in aluminum and tungsten. © All of the range distribution curves 


exhibited "tails" of deep penetration. This phenomenon was attributed to 


channeling of ions in preferred directions in the crystal. Davies, et. al. 





stated, 
The magnitude... (of crystal lattice effects)...particularly 


emphasizes the need for using either single crystals or 
amorphous materials in future range experiments. 8 


9 


In a more recent study the same group~ studied single crystals of 


tungsten using 1 to 20 kev Xel29 ions, Tungsten has a body center- 
ed cubic (b.c.c.) structure in contrast to aluminum and Cu which have 
face centered cubic (f.c.c.) crystals. General results of the single 
crystal studies were: 
...crystallographic effects are even larger than those found 
previously in Al; furthermore the results are in qualitiative 
agreement with the theoretical prediction for b.c.c. lattice, 
i.e. the (100) and (111) are the most favoured directions for 
Hianneling and the (110) and (112) are less favoured. This 
contrasts with the f.c.c. structure, where both theory and 


experiment find (110) the most favoured and (111) one of the 
less favoured.? 


Lutz and Sizmann - 0 


bombarded copper (f.c.c) single crystals with 
Kr ions of 10 to 140 kev. energy. Their results agreed in general with 
Davies‘ single crystal data, since they showed that 110 channel shots 
had deepest penetration. Further, the range distribution curves showed 
the characteristic tail of deep penetration (up to 4000A) for 110) i 
and 100 channels. Their 25 kev and 140 kev range distribution curves 
(as well as some of Davies' curves) are shown in figure 1. The experi- 
mental procedure used by Lutz, et. al. 10 employed the removal of thin 


layers of bombarded crystal by controlled sputtering with non-radio- 


4 





active krypton ions. Theradioactivity was then measured. 

In conjunction with the monocrystal experiments described above 
and computer simulation studies to be discussed, Lehmann and Leib- 
fried! have studied the behavior of primaries in channels analytical- 
ly. Their work used copper primaries in copper, as have the majori- 
ty of computer programs. Their results, oriented toward long range 


penetrations, in part were: 


Two potentials are used: an exponentially screened Coulomb 
potential after Bohr, used also in the machine calculation and 
thought to give an adequate description for relatively high en- 
ergies and small interatomic distances; anda purely exponen- 
tial potential after Born-Mayer, better suited for relatively low 
energies and large atomic distances. The maximum ranges are 
very large, for 10 kev in the order of 103 lattice parameters for 
the Born-Mayer potential and up to 10’ for the Bohr potential. 
Presumably, the Born-Mayer potential is a better description for 
these events. The investigation is confined to motions near 
the channel axis. 





2. COMPUTER SIMULATION STUDIES 

This project studies the interactions of energetic atoms in crys- 
tals through the use of a digital computer simulation program. This 
section will present a general orientation of the use of computer sim- 
ulation programs and will discuss the current work in this area. (A 
more detailed discussion of the program and model used in this study 
will be presented in section four and in Appendix II.) 

Simulation programs for any system have the following three es- 
sential ingredients. (1) The system or phenomenon has a status. 

This status is described in terms of attributes. In the case of a crys- 
tal, such attributes as position, velocity, mass and energy are essen- 
tial. (2) The status is changed or modified by an event. Events cause 
the atoms to leave their sites, gain or lose energy, etc. (3) The 

third essential element of a simulation program is time. The status and 
events are strongly intertwined in a time continuum. The control of the 
occurance of events is critical. 

Three basic approaches to simulation programs of atom-crystal in- 
teractions that have been used to date are: (1)ramdom target atoms(i.e. 
lattice effects are neglected), (2) binary collision ina lattice, and 
(3) many body collision ina lattice. The binary collision assumes in- 
teraction with one atom at a time, while the many body approach (called 
n-body) simultaneously accounts for forces from several crystal atoms. 


One of the earliest computer studies was made by Gibson, et. al. - 





Their model placed the copper atoms on sites in the f.c.c. crystal 
and used an n-body approach. To simulate an infinite crystal, add- 
itional forces were added to the surface atoms. Because the poten- 
tial functions and their corresponding force functions were not 
known precisely, Gibson, et.al. used "a simple central difference 
procedure...which gives reasonable accuracy. "12 Their techniques 
used f=ma in an ;*iteration procedure. This study showed evidence 
of chains in the 100 and 110 directions. A chain is defined as the 
propagation of energy along close packed rows of atoms ina crys- 
tal without mass transport. Interatomic potential, surface effects, 
lattice defects, and other low energy (<400 ev) phenomena were 
studied. 

The random lattice mentioned above was used inan earlier study 
by Oen, et. al. for copper primaries in solid copper. Their pro- 
gram assumed independent binary collision (accomplished by an im- 
pact parameter restriction of one half the nearest neighbor's distance), 
classical scattering, and the Bohr hard sphere approximation to the 
Bohr potential. They concluded: 

It is found that neither the hard sphere approximation nor the 

inverse r squared approximation to the Bohr potential is partic- 

ularly good. To obtain correspondence with experimental re- 
sults it is found that the Bohr screening length must be in- 


creased as the atomic number of the interacting atoms increas- 
es.l 





A more recent program by Robinson and Oen!8 (hereafter referred 
to as RO) used the binary approach again, but replaced the random 
distribution with a copper f.c.c. lattice. This program studied sev- 
eral potentials: the Bohr, eroded Born-Mayer, and truncated Born- 
Mayer. The Born-Mayer potential selected by RO had the same para- 
meters as the Gibson II potential, but incorporated four different 
truncation values. 

Significant among their conclusions were the following: 

(1) Ranges of primaries were strongly dependent on direction, 
and in f.c.c. crystals the order of range was (011)>(001)>(111)x 
isotropic. 

(2) Larger range values were due to many glancing collisions that 
confined the primaries to regions of low potential. This phenomenon of 
deep penetration is called channeling. 

(3) The combination of Bohr potential and 1000 ev primaries im- 
pacting in the (011) direction produced several primaries that moved 
from channel to channel. That is, RO found that the primaries showed 
large spread as they penetrated the crystal. Such spread was attributed 
to the nature of the Bohr potential, G6ee figure 2, which is a projection 
of their (011) channel trajectories onto the (011) plane. Note that three 


of the four tracks show large amounts of spread.) 





A more recent n-body model was used by Gay and Harrison!4 


(hereafter referred to as GH). They used an iteration procedure, 

but one that has important differences from Gibson's approach 

(see Section 4). Further, the crystal size was smaller than Gibson's, 
and the program was simplified to reduce computer time. A compari- 
son between binary collision characteristics and n-body collision 
approximations was performed. In general terms, the report showed 
that binary collision approximations approach true n-body collisions 
for small impact parameters. But at larger impact parameters the 
binary collisions have less energy transfer from primary to target, re- 
duce the primary scattering angle and slightly increase the recoil 
angle. Their work provides specific quantitative results for interac- 
tions in terms of energy exchanged, scattering angles, and recoil 
angles. 


A modified form of the GH model was used in this study. 





3; STUDY OBJECTIVES 

This study of energetic atoms in crystals was performed using an 
n-body model approach, incorporated in a computer simulation pro- 
gram. Specifically the study objectives were: 

(1) The verification and more explicit delineation of the spread 
phenomenon found by RO, ue (figure 2). As mentioned earlier, three 
of the four copper primaries shown in the figure drifted considerably 
in the lateral direction as they penetrated the crystal. 

(2) The determination of the range distribution curves for copper 
atoms penetrating the (011) surface of f.c.c. copper crystals. Spec- 
ific attention was placed on the effect of impact point on primary 


range and spread. 


10 





4. SIMUIATION MODEL 
A. Model Discussion 

The original purpose of the GH program was a comparative study of 
binary versus n-body collision processes. Their model used a 63 atom 
f.c.c. microcrystallite. These atoms were capable of movement but 
could not dissipate energy except by further collision. No surface ef- 
fects were allowed ; which was clearly justified by the short time span 
of the program. An outer shell of immovable atoms was used to simu- 
late a continuation of the lattice and test for complete event contain- 
ment. (This shell is eliminated in the current program), 

The objectives of the GH study were oriented toward analysis of the 
primary atom and one "target" atom. Concentration on the primary's 
history allowed a streamlining of the program that removed unneeded 
operations and reduced program run time. Essentially, the computation 
and recording of the motion of the remaining crystal atoms were deleted. 
This deletion was possible since the path of a primary is rarely exposed 
to feedback from target interactions. Inanalogy, if the path of radia- 
tion damage is represented by an expanding "tree of interactions", then 
the primary is usually spearheading the top branches. 

Another result of the concentrated study ona primary anda single 
target was the reduction in crystal size. It must be made clear that the 
reduction in crystal size has been carefully studied in earlier work by 


Harrison. a Considering the uncertainties of the crystal potentials and 
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the use of iteration techniques, the neglect of more distant crystal 
atoms is not significant in regard to primary performance. The crystal 
reduction was carried out experimentally through trial runs on the com- 
puter. The size was reduced until the program run time was optimized, 
but iteraction results still correlated closely. 

The unit crystal used in this study regenerates itself in one of two 
ways. The first form constructs the crystal in the y direction only 
(figure 3) while the second, slower program can construct the crystal in 
any direction (figure 4). As a result of these features, the small unit 
crystal is repeatedly remade in front of the primary. The fast running 
program is used in those cases where the primary is not expected to ex- 
perience severe collisions nor deviate greatly from its original track. 
This type of program is called Chan2. The second program is used in 
areas when the primary does experience large amounts of spread and is 
designated Chan 3. Figures 3 and 4 indicate the unit crystal sizes and 
their orientation. All primaries hit the 011 surface orthogonally (i.e. in 
the y direction). The crystal is positioned by the program with 011 sur- 
face in the x-z plane and with the positive y direction coinciding with 
the primary track at impact. A majority of the runs were made ina tri- 
angular impact area. Several check runs were made in the adjecent triangle 
and a few elsewhere. Detailed discussion of the impact triangles is given 
in the next section. 


° 


The basic building block of a copper crystal is a cube measuring 3.614A 


ona side. Tne simulation program uses units of length of 1.807A called 


be 





lattice units. 





The program is written in Fortran and Symbolic Fortran (an assembly 
language closely related to machine language). All runs were conducted 
on a Control Data Corporation 1604 computer which has a storage capacity 
of over 32,000 words and six index registers. A CDC 160 computer was 
used to automatically plot graphic projections of the primary trajectories. 
Refer to figure 5 for an example of such tracks projected onto the 011 plane. 

The program constructs the lattice (i.e. assigns positions to each atom), 
commences the run, at designated time intervals computes the forces on the 
primary and secondaries, adjusts their motion accordingly, prints required 
information at set intervals, and "shuts down" when any of three conditions 
exist. The program terminates when (1) the energy of the primary goes be- 
low 25ev, (2) a preset time limit is exceeded (measured in program cycles; 
10 cycles ~w time for the primary to travel one lattice unit), or (3) the primary 
leaves the lattice een for Chan 2 programs only). The determination 


of this cycle factor received considerable study by Gay? 


and was optimized 
experimentally. 

As may have been evident in the brief program discription above, the crit- 
ical phase of the program is the method of computation of the primary pathas 
it experiences crystal forces. In general, the "events" of this simulation 
program are the primary crystal interactions. 

The interaction routine used by GH is still used in this model. It consists 


of a double iteration procedure that calculates the average force on an atom 


over a small distance of travel. This force is used to advance the atom. 
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meticon. states; 


The unbalanced force ...(on the primary)...is an average force 
calculated by a double iteration procedure as follows: (1) assume 
an atom at position 1 with velocity 1; (2) calculate the total force 
on the atom as a result of all the other atoms in the lattice (this 
means normally only about 8 ~ 10 nearest atoms...); (3) call this 
calculated force, force 1, and use the equation of motion to move 
the atom to a temporary position, position 2; (4) now repeat the 
force calculations for position 2, call this force 2; (5) go back to 
positional, and use the average of force 1 and force 2 to move 
the atom to a new position, position 3. Procedures 1 through 5 
constitute one "time step". 14 


The equation of motion used in the iteration procedure differed in two 


12 


ways from the equation used by Gibson” et. al. Although derived ina 


different manner, the GH equation can be compared more easily to the 
Gibson form by starting with the latter and indicating the modifications. 
The Gibson form of the equation (in the x direction only) is: 

X; (t + At) = X; (t) + Atl vy (t+ At/2)+ m7!. F(t). At| 
where X; = initial X coordinate, t = time, At = time increment for change 
of position, vj = initial velocity, m= mass, and f= the initial force. 
Because At is a small number compared to t the term Vj (t + At/2) may be 
expanded ina Maclaurin series. Neglecting all terms after the first , the 
following form is obtained: Xj (t+ At) = xj + At iV; (t) + At/2. F(2)/M| 
In the GH derivation, the force term is an average force where 


ie a Fit lhe 
is 5 
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and the following expression results: - 

Xj (t+ At) =Xji+ At | Vi ay + St EW | 
Thus, the GH expression resembles the Gibson expression but utilizes 
a more realistic force term. 
B. Potential Functions 

Two factors dominate all other considerations in the simulation of pri- 
mary interaction in crystals. The first is the general type of model attempt- 
ed. Specifically, does the program use binary or n-body collisions? Isa 
crystal or a random atom arrangement used? The second factor is the se-- 
‘estion of a potential function. Many functions have been proposed; none 
have been singled out as best. Until good correlation between experiment 
and computer simulation is made, the proper potential function will remain 
uncertain. Ina sense this report and all such simulation studies are an at- 
tempt to locate the correct function for a particular substance. 

Two basic types of potential function were used in this study: Born- 
Mayer and Bohr. The Born-Mayer function used has three parameter varia- 


12 


tions introduced by Gibson et. al. and shown below: 


@¢=Aexp |. p (r-r,) /x, | ; 
where r = distance in lattice units, 


r, = nearest neighbors distance (2.555A for copper), 


°o 
1, = lattice unit (1.804A for copper), 


po = input parameter, and 


> 
l 


input parameter (ev). 
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The form shown above is convenient for computer calculation. The 


three sets of parameters selected by Gibson are: 


Potential A(ev) 0 
1 . 0382 1Lono7 
II 951 13.00 
IIT - 1004 10.34 


To speed calculation, the Born-Mayer force functions were truncated at f=1 oe 


newtons (corresponds to approximately 3A). 
In similar format the Bohr potential function is: 
ee | exp (tc, ay-r/ay), 
where r = distance in lattice units, 
a,= input parameter in lattice units, and 
C,= input parameter in ev. 
The two parameters are ; cy = aWos ioe and a,= .0674lu. These para- 
meters are the ones used in the RO report. Further, the corresponding force 
truncation value was used so that the Bohr potential functions and force func- 
tions in this study are identical to Robinson's. (The truncation of the force 
function occurs at f = 3.9x107” newtons). 
The lattice has zero kinetic and potential energy prior to primary impact. 
In other words no vibrational energy is present. 
C. Binary vs n-body model 
This study made use of an n-body simulation program in contrast to the 


ee ogy 


simplified binary collision approach of Robinson, et. al. The in- 
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herent danger of superimposing independent collisions is that simultaneous 
collisions will not be properly accounted for. For purposes of fast comput- 
er runs, the binary program has the greatest advantage, but its justifica- 
tion requires better experimental verification than is yet available. The 
n-body program is inherently more accurate because of its more thorough in- 
teraction routine. Given proper parameters, it should yield the more accurate 


data. 


I, 





9. EXPERIMENT PROCEDURE 

Computer simulation studies are experiments in the accepted sense of 
the word. They are trials made to confirm or disprove some suggested 
phenomenon, and are subject to most of the usual experimental inaccuracies. 
Systematic errors occur because of improper model simulation or necessary 
oversimplifications of the model. Random errors caused by computer mal- 
function should be rare. Provided the program and its input data are not 
modified, the computer should reproduce results perfectly. (Refer to subsec- 
mon C.). 

Simulation studies are sensitive to parameter variation to the extent that 
some computations can become meaningless if not properly controlled. For 
example, the combination of high energy primaries, Bohr potential, and sel- 
ected impact points produce an interaction so weak that the program cannot 
maintain the energy check. As a result the primary could penetrate infinite- 
ly. Such extreme situations have been avoided and are not included in the 
data of the report, but they do point out the danger of blind faith in computer 
"answers", 

The remainder of this section will present the experimental approach used. 
A. Determination of impact area and sample size 

Symmetry considerations reveal a ne impact rectangle shown in figure 
3 or 4 represents all possible combinations of impact possibilities. That is, 
the rectangle is a "unit area" that, if repeated, will completely describe the 


011 surface. Further, with one exception, the rectangle can be split down 
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the diagonal into two triangles that have point correspondence in regard to 
primary interaction. Consider impacts at pointsA’and C’{the 90 degree an- 
gles). As indicated in figure 4, atom row AB is nearer the surface than row 
DC and therefore "sees" the primary sooner. Now, a primary, or "bullet", 
striking the crystal at A’ "sees" these top plane or near atoms (row AB) at 
a distance of 1.28A . However, a bullet hitting at point C’ sees the top 
plane atoms (row AB) at a distance of 1.807A . Obviously, this difference 
effects the first interaction and subsequently alters the bullet's path and 
future collisions. It will now be shown that this difference in the two 
triangles, will cause a mild perturbation in range and spread data, and is 
not considered serious for purposes of this study. (The triangular impact 
area is desirable since it reduces the sample size necessary by one half). 
A graphical correlation analysis of two critical parameters (range and 
spread) was performed to insure the validity of the impact triangle. In 
each triangle 14 points were chosen and the range and maximum spread 
values recorded. All points were obtained using 1000ev primaries and a 
Gibson II potential. Since this combination has a strong primary<-crystal 
interaction, discrepancies between the two triangles ehoula be maximum 
(figure 7). In general, good correlation was revealed. Note that the range 
correlation values show less correlation for deeper penetrations. This is 
expected since those primaries that impact further from the “corner atoms" 
have smaller energy loss rates and more chance to deviate. The spread 


terms have less correlation at small values. However, these terms are 
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Sufficiently small so that .2 or .3 lattice unit variation between triangles 
does not effect the spread contour diagrams (to be discussed later). Asa 
result of these considerations, triangle ABD was selected as the primary 
impact area. The error introduced by this assumption is small and is not 
known precisely. In future studies one area of further investigation should 
undoubtedly be a more thorough check of the rectangular impact area and its 
effect on results. 

From the grid shown in figure 8, 14 impact points were selected and 
are labeled one through 14. For some of the runs these original 14 points 
were supplemented by 12 more. These points are labeled in the figure as 
1A through 12A. The supplemental 12 point array was chosen to verify the 
selection of the basic 14 point sample size. Both the Gibson I[J-1000 and 
59000ev studies were initially made on 14 point grids and later supplemented 
by 12 more data points. The range contour surfaces remained regular when 
the new points were added, and showed no anomalies (figure 9). Further, 
the range distribution curves for the 14 point and 26 point grids correlated 
very closely (figure 6), and established the soundness of 14 point grids. 
These range curves were constructed on an area rather than a point basis. 
That is, cross sectional areas of the contour surfaces (figures 9 through 12) 
were measured at specified depths of penetration. The ratio of this area 
to the total triangle area gives the percent of primaries stopped. In figure 
9 the shaded area represents those primaries stopped within 100 lattice units. 


The ratio of this area to the total area provides the percent stopped at 100 
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lattice units. (The range curves were constructed from larger contour 
diagrams than those shown). The area method yields more accurate pene- 
tration curves than woulda point system. Overall justification of a 14 
point sample is simply the general regularity of the results. 

B. Impact points A,B,C and D. 

In addition to the 26 point grid, four other points were chosen on the 
triangle. These points coincide with the impact points in figure 3 of the 
RO study (three of the four points showed spread). These points are 
shown in figure 8 of this report and are labeled as A,B,C and D. 

As an initial step in this study, the parameters of the RO "A,B,C 
and D shots" were duplicated as closely as possible. Later, several 
other combinations of bullet energy and potential were used. Table 1 
shows the number of grid points investigated under each combination of 
parameters. As will be shown in the next subsection, many of these 


points were rechecked with additional runs. 
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TABLE 1 








Impact Points Studied By Parameters 


ro! Points l4 Point 12 Point Symmetrical 
AG BEC 7) Grid Grid Triangle 


Energy Energy Energy Energy 
Potential 1000 5000 1000 5000 1000 5000 1000 


Bohr 4 Point C 14 is 8 6 - 
4 Point © 14 14 12 12 14 
4 = 14 - = - - 


Gib II 
C. Reproduceability of runs. 
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As stated earlier reproduceability problems should not occur in com- 
puter studies. Unfortunately, the program - computer combination used 
in this study does fail to reproduce periodically. Asa result a large 
number of reruns were necessary. The sample was random, but was 
later checked to insure that most combinations of energy, potential and 
run length were represented. The sample size was approximately 15 per- 
cent. It was found that the probability that a run would be reproduced 
exactly was 85 percent. Further, once a run was confirmed by a rerun, 
the probability d its correctness could be assumed to be one. The lack of 
perfect reproduceability has not been determined. However, it was found 
that those runs that deviated seemed to get out of phase, but did not vary 
greatly in amplitude or distance traveled. Generally, the variation of am- 


plitude (i.e. spread) and depth of penetration were less thana few per- 
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cent. In one case the depth of penetration varied by about 50 percent, 
but the depth of penetration in this case was very short ( <8Lu). 
D. Program output. 

The output will be briefly discussed here, so that subsequent dis- 
cussions of data will have more meaning. The program as used in this 
study was designed specifically to examine only the primaries' per- 
formance. Output is made every ten cycles. Ten cycles is the time 
necessary for the bullet to travel approximately one lattice unit (1.807A). 
Each output yields the following data about the primary: 

(1) x,y, and z distances (measured in z lattice units) from the im- 
pact point (a "z lattice unit" is 1.414 times greater than either the x or 
y lattice units), 

(2) X,Y, and z velocity values in "z lattice units" per second, 

(3) kinetic energy, potential energy, and total energy in electron 
volts, and 

(4) ‘real elapsed time in seconds. 

Except for a few 6000 cycle runs all of the simulations used a 2000 


cycle cut off. 
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6. RESULTS AND CONCLUSIONS 

The data collected from the printout sheets discussed in section 
five were condensed and tabulated as follows: (1) primary final energy 
(either approximately 25ev or the energy at the 2000 or 6000 cycle cut 
off), (2) depth of penetration in lattice units, and (3) the maximum 
spread. Recall that spread is the perpendicular radial distance (in 
lattice units) that a primary deviates from its initial straight line tra- 
jectory. 

The type of tracks reported by RO (those of large penetration and 
spread) suggest the need for a word descriptive of their nature. There- 
fore, in this report, the term wander will apply to all primaries that ex- 
ceed two lattice units of spread, but whose energy loss rate does not 
exceed ee 

Four combinations of potential and energy were studied: Bohr and 
1000ev primaries, Bohr and 5000ev primaries, Gibson II and 1000ev 
primaries, and Gibson II and 5000ev primaries. Since the data was 
collected from the impact point grids, it is best displayed using a pic- 
torial representation of the grid. The penetration values are shown at 
their respective grid impact points and can be visualized as a penetra- 
tion contour surface (figures 9 through 12). The pertinent data from the 
spread contours was summarized and is discussed below. 


A. Spread results 


The spread values for each impact triangle are a function of primary 
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energy, potential and impact point. Separate spread contours are in- 
cluded in the report, and are located on the same page as their corres- 
ponding penetration contours. 

The results of the spread data collected are summarized below: 

(1) Wander, as observed in the Bohr 1000ev primary case by RO, was 
found to exist for several other combinations of primary energy and poten- 


tial function. 
TABLE 2 


Percent of primaries that exceeded 2LU spread 
and Percent of primaries that wandered 


POTENTIAL 


Energy nergy 


1000ev 9000ev 1000ev 5000 ev 


Percent 43 
>2LU Spread 





The result of RO can be generalized: When parameters are such that 
penetrations vary from deep to short values over the surface of impact 
triangle, then someof the primaries will wander. Note that the Gibson 
IIT - 1000ev series of runs does not satisy this criteria (refer to table 2), 


since all the primaries had short penetrations. 


From another viewpoint, if the strength of interaction between primary 
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and crystal varies over the impact grid from hardto soft collision, then 
some primaries will wander. Thus, the presence of wander is not simply 
a result of the Bohr potential, but is a function of the strength of inter= 
action. 

(2) A study of the spread data (figure 9 through 12) reveals that the 
primaries that wander are located in a particular region of the impact 
triangle. Note that the primaries impacting in zone A, figure 14, have 
large spread values while those in zone C have large penetrations. It 
is not surprising that the primaries that wandered were located at the 
intersection of zones Aand C. The overlap of the two zones varies some- 
what with the parameters of the run and cannot be firmly located. How- 
ever, the important result is that wander is a function of impact position 
and is simply the transition region from the short range penetration zone 
to the deep penetration impact zone. 

(3) The outlined rectangular area in zone C (figure 14) has primaries 
that have low spread. Any atom impacting this area will have large pene- a 
tration and will not exceed.8Lu spread. The primaries in the remaining 
area of the impact triangle have at least lLu of spread except for a 
small area near point 11; the midpoint between the atom crystals. This 
midpoint is a symmetrical force region for the primaries and little deflec- 
tion or spread takes place. Further, in the Bohr potential cases, the force 


curve is near cut off so that the forces on the primary are negligible. Two 
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regions of large spread occur between the .3 and .5Lu radii of the cor- 
ner atoms. These are the zones where the primaries wander. 

(4) Primaries 2, 5A, and 8A in figure 15 wandered. As shown by 
the "history curves" of the runs, the rate at which energy is lost is us- 
ually uniform until approximately 1000-1200ev and then the energy loss 
is more rapid. The general appearance of the family of curves is uniform, 
and one of gradual transition from short to long penetration as the impact 
points move from the corner atoms of the triangle to the channel area 
near the right angle. The short penetration curves display less regularity. 

The gap in the region of 180A to 300A penetration is not an amonaly. 
If primaries were concentrated in the region of rapid penetration change 
(i.e. the "hills" of the penetration contour) then this gap would be fill- 
ed. Though the crystal effecton range and spread seems logical, the 
rigid consistency and smooth transition is surprising. For example, if 
more impact points are studied, the occurance of occasional penetrations 
in the range of 300A to 600A that do not wander can be expected. 

In brief summary, primaries do wander if the impact point is located 
in a transition region from hard interaction to soft interaction. The "wan- 
der phenomenon" is part of a smooth change from one type of interaction 
to another. 


(5) The four impact points of the RO study that wandered were investi- 


gated (figure 2). Their parameters were duplicated (Bohr truncated poten- 
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tial, 1000ev) to insure a valid comparison. In addition, the Gibson II 
and Gibson III potentials were used. The selection of the exact para- 
meters in the Bohr case provided a comparison between the RO binary 
model and the n~-body model. 

Of the four points, only point C displayed wander and it did so with 
the following combinations of paramenters: (1) Bohr and 1000ev pri- 
maries, (2) Bohr and 5000ev primaries, and (3) Gibson and 5000ev pri- 
maries. As seen in figure 14, points A, B, and Dare in zone C and 
are typical of that zone Since they displayed deep penetration and less 
than .8Lu spread. The probable reason for the disagreement between 
the RO study and this study is the use of the n-body model. (Refer to 
the next subsection for detailed comparison of the results obtained from 
the two models.) 

At present, none of these results can be verified experimentally. If 
the potential functions studied are not reasonable, the conclusions dis- 
cussed above could be radically change. As will be seen in the next sec- 
tion, the potential functions selected (particularly the Bohr) may be too 
weak. If so, the amount of "wander" found in future studies using "hard- 
er" potentials may decrease or disappear. 

B. Penetration results 

The penetration results are shown in contour diagrams (figure 9 through 

12) and the integral penetration curves of figure 6. Penetration contours 


for all runs are regular and well defined. The following points are reveal- 
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ed from the data: 

(1) Impact points 4, 4A, 7 and 7A, which are located in the 110 
channel of the crystal have maximum penetration. The Gibson II- 
1000ev contour is complete and reveals the smallest and lowest 
peak of maximum penetration. The Gibson I[-5000ev contour is virtual- 
ly complete (four runs truncated at 600Lu) and is shaped similarly to 
the Gibson II 1000ev contours. However, the contour height is much 
greater and the region of very low penetration (below 20Lu) is much 
smaller. 

Both Bohr contours were truncated at 2000 cycles. In comparison to 
the Gibson results, the areas of low penetration are smaller and the rate 
of change to deeper penetration more rapid. Excessive computer run 
time curtailed the construction of complete contour surfaces for the Bohr 
potentials. 

The Bohr-1000ev contour and the Gibson II-5000ev contour are simi- 
lar in size and shape, which indicates that these combinations of ener- 
gy and potential produce similar interactions. 

(2) As the energy of the primary is increased the depth of penetra- 
tion increases; also, the amount of penetration varies with the potential 
function used. 

The n-body integral penetration curve (5000ev and Gibson II poten- 


tial has much larger penetrations than the RO results (figure 6). Un- 
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fortunately, none of these simulation curves can be easily compared 
to the experimental data of figure 1, because of the large difference 
in the substances studied. Important conclusions are available from 
the results, however. 

The 5000ev curves in figure 6 of RO and this study agree fairly 
close in the region of small penetration (up to 120A). Thus, it is pro- 
able that the n-body program agrees closely with the binary model for 
severe interactions (i.e. small impact parameters below .3Lu or 5A). 
A thorough study by GH established the close agreement between single 
binary collisions and n-body collisions at short impact parameters. 
Their study investigated a wide range of energies and several poten- 
tials (including the Gibson II). 

(3) In the region of deep penetration the two models disagree. 
Nevertheless, an important fact was revealed when the simulation pene- 
tration data were compared with the corresponding potential curves. 
Refer to figure 13 which shows the Gibson II potential discussed here, 
the Bohr potential, and a new potential to be explained in the next sub- 
section. As seen in the diagram, the three curves in question use the 
basic Gibson II potential function. However, RO curves B and C are 
truncated at much larger potential values than the potential of this study. 


As seen from the integral penetration curve and these potentials, as 


the distance to the truncation point is increased, the depth of penetra- 


30 





tion of the primaries increase . (This effect is opposite of what might 
be first supposed). The most plausible explanation probably lies in 
the subtle effect of the small potential curve slope (hence small force 
values) that occured in the region beyond the truncation of the RO 
curves. These small forces confine the primaries to the more open 
crystal regions. As an analogy, the deeper penetrations can be com- 
pared to a rock skipping over the surface of water. It makes many 
weak impacts that confine its movement. This is the case of the n- 
body potential. The more radically truncated potentials of the RO study 
probably cause abrupt interaction when a primary moves from a region 
of no potential to one of large finite value. Thus, deep penetration of 
primaries still occurs in the RO model but toa much smaller extent. 

In summary, one significant factor causing the deeper primary pene- 
trations of this study was the use of a potential that was truncated at 
large radii. The use of the n=body model also contributed to large pene- 
trations and will be discussed shortly. 

(4) As shown in figure 13, the RO study also used an eroded poten- 
tial that decreased to zero at 1. Sah The corresponding integral pene= 
tration curve is similar in appearance to the n-body penetration curve of 
this study, but has slightly lower penetration values. Thus, the poten- 
tial function if eroded or truncated at large distances tends to produce 


Similar large penetration depths. In contrast, potentials that have signifi- 
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cant values at truncation reduce penetration depth. 

(5) Another factor, which undoubtedly caused the deeper penetra- 
tions found in this study, is the n-body model. It’s use insures a 
continual superposition of all significant crystal forces; and consequent- 
ly, the net force on the primary is usually smaller and produces less 
scatter. Thus, the wander detected by RO for impact primaries A,B, and 
D, but which was not verified using the n-body model, may result from 
the insufficient superposition of forces in the binary approximation collis- 
ion process. This helps clearify the results in the previous subsection 
concerning wander. 

C. Recommendations 

Two proposed simulation studies are presented in this subsection. 
Both are attempts to obtain verification with experimental data. Recent 
studies by Peereinicena provide potential functions for tungsten and the 
inert gases. This development, in conjunction with the experimental 
work of Davies, provides, for the first time, the possibility of positive 
experimental verification of computer simulation programs. 

It is recommended that the n-body program incorporate the new poten- 
tial function developed by Abrahamson, and that it be used to simulate 
the Davies experiment of “eno in tungsten. 

Although the experimental results for copper primaries in copper are 


not available, it is felt that the simulation penetration values are too 


32 





large (especially the n-body) and should agree more closely with the 
25 kev - ica in Cu experiments and the cna in Cu experiments at 
Skev. In this regard RO have used "Cu" primaries of triple mass 
ole which did reduce penetration ranges (fig. 6), but not suffi- 
ciently to agree with the results of 25kev Geo in Cue 

In spite of the probability that the n-body integral penetration curve, 
obtained in this study, is further from experimental fact than the RO 
curve, the n-body model is an inherently more accurate approach. This 
conclusion is based on previous discussion which reduces to the essen- 
tial fact that the n-body model is more realistic and, hence, a better 
simulator, Equally critical, the results of the RO study were obtained 
using severe truncation values (fig 13). Truncation of the force func- 
tion at significant values is unrealistic and probably leads to incorrect 
results. To obtain closer verification of experimental results, stronger 
potential functions may be necessary in future n-body program studies. 
(It is possible that experimental integral penetration curves could be 
verified by most models by sufficient adjustment of the potential 
functions. Thus, experimental verification ina single area will probab- 
ly not prove sufficient, in the long run, to establish a particular poten- 
tial function). 

Preliminary calculations indicate that the potential function of 


Abrahamson (Xe in W) is significantly stronger than the Gibson II poten- 
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tial. If so, it is likely that this potential will produce excessive 
duction in primary penetrations and its use ina simmulation program will 
not yield close correlation to Davies experimental work. 

It is also recommended that the verification of Davies experiment 
be attempted using a "compound" force function. This function should 
follow the Gibson II force curve at small impact parameters where 
good agreement is obtained at present. At approximately 1. ip the 
curve would change toa different exponential form that exerts much 
larger long range forces on the primaries. This function will reduce 


deep penetrations and agree more closely with experimental data. The 


parameters of Davies Xe in W should be duplicated for this study also. 
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Fig. 9. CONTOUR CURVES. Gibson II potential and 5000ev primaries. 
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Fig. 12. CONTOUR CURVES. Bohr potential and 1000ev primaries. 
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APPENDXK II 
PROGRAM DISCUSSION 

A. Program Explanation 

The descriptions given below are to be used with Appendix I (Pro- 
gram Listing), and are specifically applicable to programs Chan3, Bull- 
et3, or G Chan3. 
BOX 1 

Listing of dimension, common, and Format statements. Provide graph 
plot data needed such as graph size, scale, and titles. Read force and 
potential function parameters. Read primary atom parameters. Evaluate 
miscellaneous constants. 
BOX 2 

Designation of nearest neighbor lattice distances. Construction of 
microcrystallite. Three nested "DO loops" that extend to statement 63 
are used to locate and number the atoms in the crystal. 
BOX 3 

Initialize the velocity and the position of the primary. 
BOX 4 

Store initial positions of all atoms. Initialize various terms (i.e. 
set=0). Print variable quantities such as primary energy, direction 
cosines, mass, etc. which comprise the output heading. 
BOX 5 


Compute several constants given DT, CVD, PTMAS and PGMAS. 
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BOX 6 

This entire box is a "DO loop" onall crystal atoms. Set all forces 
to 0. Check coordinate distances from primary to lattice atoms and 
neglect if greater than ROE. Compute the square of the "vector" dis- 
tances between the primary and lattice atoms and neglect if greater than 
(ROE)? - Compute the force that corresponds to the vector distance and 
compare ittoFM. Neglect if less than FM. Compute forces on the 
primary and lattice atoms. 
BOX 7 

If index=1 jump to box 8. Compute temporary position of the primary 
and then the rest of the crystal atoms. All atoms are advanced at their 
initial velocity for DT seconds, while experiencing forces found in box 6. 
BOX 8 

(Return to box 6 and recompute the forces on all atoms. These values 
are added to the initial forces yielding force terms twice the average. 
Since index=1 on reaching box 7 a jump is made to statement 260 in box 8). 
The averages of the initial and temporary position forces are used to compute 
the final positions of the primary and lattice atoms. (See page 14 of text). 
BOX 9 

Perform tests to see if program output is necessary or if the program 
should rebuild the crystal. If the latter is true jump to the appropriate re- 


generation box. 
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BOX 10 

Calculate the distance from the primary to each lattice atom; if less 
than ROE calculate the potential energy of the interaction and store. 
BOX 11 

Calculate the quantities needed in the output statement (i.e. pri- 
mary position, kinetic energy, total energy, and real elapsed time). 
Write the output statement. Check to insure the potential energy is less 
than "Bitty". Terminate program if the potential energy is less than "Bitty" 
for three passes. 
BOX 12 

Check current program cycle. If at "cut off" cycle, terminate the 
program. 
BOX 13 

This section regenerates the microcrystallite in the+y direction. 
BOX 14A 

Regenerate microcrystallite in the -X direction. 
BOX 14B 

Regenerate microcrystallite in the +X direction. 
BOX 15A 

Regenerate microcrystallite in the -Z direction. 
BOX 15B 


Regenerates microcrystallite in the +Z direction. 
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BOX 16 


Call the graph plot subroutine (called Draw). The following parameters 


are necessary: 


Variables 


II 


GRX 


GRZ 


C 


C 


LABEL 


ITITLE 


EXSCALE 


YSCALE 


0 


15 


ERROR 
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Explanation 


Number of points to be plotted 
X coordinate 


Y coordinate (simulation program 
Z coordinate) 


Number of curves per graph 
Point or curve plot 

Type End by the end of curve 
Graph title 

Graph scale in X direction 
Graph scale in Y direction 


Distance of x axis from bottom 
of graph 


Distance of y axis from left edge 
Mode of x axis (automatic) 

Mode of y axis (automatic) 
Graph width 

Graph heighth 

Place grid on graph 


Indicates if previous plot was OK 





B. Definition of Variables: 





Term 
Bitty 
BE 
ba, BY, BZ 


COX, COY, COZ 
CVD=1.8075x107!9 


Cve=1.60x10 29 


CVM=1.672x10 2" 


DIST 


DT 
DTI 
DTOD 


DRX, RRY, DRZ 
DXS, DYS, DZS 


Bat, OYT, DAT 
El=25.0 


EXA, EXB, EXC, EX1 


EV 


Usage 


a small increment that is compared to 
primary potential energy (PTE). Used 
to terminate program when potential 
energy goes below increment value. 
primary energy (in joules) 


primary initial coordinate (in lattice 
units) 


direction cosines of primary 
converts lattice units to A (Cu only) 
converts ev to joules 


converts AMU to kg. 


distance between two atoms or the 
square of the distance (lattice units) 


basic time step length (seconds) 
variable scale factor for DT 
DT/CVD (a “units Converter”) 


coordinate distance of primary to 
lattice atoms (in lattice units) 


locates current microcrystallite with 
respect to original origin (lattice units) 


test constants used to rebuild crystallite 
cut off energy (in electron volts) 


input potential function parameters 
(see Gay!® thesis) 


primary energy (in electron volts) 
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Pay re, FC 


FM 


FOD 


FORCE 


PXA, FXC 


feeeeey, FZ 


GMAS 


GRX, GRZ 


Hex, IDY, (DZ 
Mell, KT 
i, ly, 12 


IXP, IYP, IZP 


eT 
HDTOM 
HDTOMB 


HGMAS 


HTMAS 


=2 


x,y, and z force components be- 
tween two atoms (in newtons) 


a small force increment used to eli- 
minate weak force terms 


force/distance (in newtons/lattice units) 


eroded force between two atoms 
(in newtons) 


input force function parameters 
(see Gay thesis 


x,y,z coordinate force terms, respectively 
(in newtons) 


primary mass (in AMU) 


coordinates of primary that are used by 
graph plot routine (in lattice units) 


constants used in lattice reconstruction 
test points in lattice generator 
unite lattice dimensions (in lattice units) 


specify volume of microcrystallite (in 
lattice units) 


used in crystal regeneration 
DT/(2*PTMAS) 
DT/(2*PGMAS) 


1/2 primary mass. units: (ev*kg)/ 
(m2/sec2) 


1/2 target mass. units: (ev*kg)/ 
(m2/sec2) 


a constant 
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MASK 


ND 
NS 
NT 


NTT 


ONE=1.0 
OPNX 
OPNZ 

PG MAS 
PKE 


POT 


Ere 


PTE 
PTMAS 
QDTOM 
QDTOMB 


ROE 


RO 
RX, RY, RZ 


BAL, RYI, R21 


used to calculate absolute value 
(removes minus signs) 


print statement increment (in cycles) 
initial print statement (in cycles) 
current program cycle 


program cycle for final print statement 
("cuts off program") 


a constant 

1 seo Gx 

IZ IesG@Z 

primary mass (in kg) 

primary kinetic energy (in electron volts) 


potential between any two atoms(in elec- 
tron volts) 


primary potential energy (in electron 
volts) 


primary total energy (in electron volts) 
target mass (in kg) 

DT/(4*PTMAS) 

DT/(4*PGMAS) 


nearest neighbors separation (in lattice 
units) 


=1/ROE 
atom coordinate (in lattice units) 


initial atom coordinates (in lattice units) 
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RXK, RYK, RAK 


BCX, oCY, 5CZ 


TIME 
TMAS 
TPO=2.1 
TPOX 
TPOZ 

VOL 

VX, VY, VZ 


ALL, ZLL 


ZE= 0,0 


temporary atom coordinates used in 
force calculation (in lattice units) 


factors to convert X,Y,Z to the same 
lattice unit basis; (values depend on 
type crystal used; i.e. BCC orf.c.c., 
etc.) 

real elapsed time (in seconds) 

target mass in AMU 

a constant 

Za aex 

Zea GZ 

velocity of primary (meters/seconds) 


atom velocity (in meters/seconds) 


test constants that determine if 
regeneration is necessary 


a constant 
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APPENDXK III 
PROGRAM OPERATING PROCEDURES 

The simulation model consists of a single basic porgram that has 
two forms; 

(1) A fast running program that can reconstruct the basic lattice in 
the y direction only (the direction of the incoming bullet). This pro- 
gram is limited to cases in which the bullet does not intereact strong- 
ly and hence doesn't leave the basic lattice in any direction except 
the y direction. This program runs in about one half the time of the 
other program and has three identical decks which are titled as follows: 
Chan, Gchan, and Bullet. For runs cut off as 2000 cycles the "Chan 
type " program requires about ten minutes. Those shots that strike 
the impact triangle near the corner atoms often terminate in very short 
time since the primary rapidly drops below the cutoff energy of 25ev. 

(2) The slower program can contain any bullet trajectory since it 
reconstructs the lattice in any direction required. The testing procedure 
to perform this function increases the program run time. This program 
has two identical decks that are titled: Chan3 and Bu_let3. (The various 
deck names serve only to expedite processing and ease administration 
of runs. The slower program (Chan3) was the third modification of the 
basic n-body program). Run times are about 20 minutes for 2000 cycle 
cases. 


The items discussed in the remainder of this section are applicable 
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to all decks: 
A. Parameter cards internal to the deck. 

Within the program deck there are three parameter cards that are often 
changed; the potential function, force function, and the force truncation 
card. The force and potential cards are the third and fourth cards in all 
decks. These fortran functions are computed froma single imput data 
card (see the next sub section). Thus the force and potential cards are 
paired to an input parameter card. Specifically, the Bohr potential and 
force cards have a single corresponding input card. The Born-Mayer 
force and potential cards become one of the three Gibson forms (I,IL or 
III) through the use of one of three corresponding input parameter cards. 
The third internal card is the force card (ex. fm = 1.0E-11) which allows 
the selection of any desired truncation force. 

B. Input data cards. 

These are the basic data cards used at the end of the deck to spec- 
ify the various parameters for a run and to set up the graph plot routine. 
The following list describes the cards in order as they would appear for 
a single run: 

(1) The first card prints any desired statement at the top of the 
print out. 

(2) This is the force and potential input parameter card and it supplies 


the follwoing values: EX1, EXB, EXC, FXC, IH. The first four quantities 
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are needed to define the function and the last (IH)prints the potential 


function in the printout. The above data cards are used once per run. 
(ChanS uses a new input parameter system). 


(3) This card provides the variable quantities that specify the 


primary parameters and general program information. The following 


1tems appear: 

GMAS = Primary mass in AMU 

TMAS = lattice atom mass in AMU 
BX = X initial coordinate for the primary 
BY = Y initial coordinate for the primary 
BZ = Z initial coordiante for the primary 
EV = Primary energy in ev 

COX = Direction cosines of primary 

COY = From the 011 plane to trajectory line 


NTT = Cut off cycle (ends program) 


NS = First print statement (in cycles) 

ND = Print increment (in. cycles) 
Recall that one cycle = .1 time for the primary to travel one lattice unit 
(1.807A). 


The next six data cards are used to title and identify the graph out=- 


put for the run. 
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(4) Potential is designated (ex. GIB. for Gibson) 
(S) Energy (ex. 5000ev) 

(6) X coordinate (ex. x = 2. 00) 

(7) Z coordinate (ex. z = 3.00) 


00) 


(8) Cosine X (ex. cox 
(9) Cosine Y (ex. coy = 1.00) 
This completes the data cards needed for single run. 

To provide additional runs cards three through nine are repeated for 
each run. An “end card" is placed last on the deck to terminate the pro- 
gram. The program is designed to run on the regular monitor operation 
service, but for runs longer than 15 minutes (the operation service limit) 
the program must be run individually. Three tapes are required, in add- 
ition to the Fortran 60 compiler tape, to run the program on the 1604 
computer. These tapes are the input, output, and graph plot tapes. 

(To perform runs without the graph plotter pull the "call graph card", 


the title cards (4 thru 9), and the card that "reads in" the title card. 


Te 
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